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FOREWORD 


This is the progress report on the research project Numerical Solutions of Three- 
Dimensional Navier-Stokes Equations for Closed-Bluff Bodies”. Within the guidelines 
of the project, special attention was directed toward research activities in the area 
of "Grid Sensitivity in Airplane Design.” The period of performance of this specific 
research was January 1, 1991 through December 31, 1991. This work was supported 
by the NASA Langley Research Center through Cooperate Agreement NCCl-68. 
The cooperate agreement was monitored by Dr. Robert E. Smith Jr. of Analysis and 
Computation Division (Computer Applications Branch), NASA Langley Research 
Center, Mail Stop 125. 
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Analysis For NACA Four-Digit Wing Sections 
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Sensitivity analysis in Computational Fluid Dynamics with emphasis on grids and 
surface parameterization is described. An interactive algebraic grid-generation tech- 
nique is employed to generate C-type grids around NACA four-digit wing sections. 
An analytical procedure is developed for calculating grid sensitivity with respect to 
design parameters of a wing section. A comparison of the sensitivity with that ob- 
tained using a finite-difference approach is made. Grid sensitivity with respect to 
grid parameters, such as grid-stretching coefficients, are also investigated. Using the 
resultant grid sensitivity, aerodynamic sensitivity is obtained using the compressible 
two-dimensional thin-layer Navier-Stokes equations. 
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NOMENCLATURE 


Q* 

x* 

R 

P 

X D 

u 

V 

xi,yi 

r 

s 

t 

R(r) 

S(s) 

I<ul<2 

k 

T 

M 

C 

C P 

Cd 


= vector of field variables 
= vector of physical coordinates 
= steady-state residual 
= vector of control parameters 
= vector of design parameters 
= horizontal interpolant 
= normal interpolant 
= physical coordinates 
= airfoil surface coordinates 
= surface grid distribution function 
= far-held boundary grid distribution function 
= grid stretching function 
= surface orthogonality function 
= far-held boundary orthogonality function 
= orthogonality parameters 
= stretching parameter 
= maximum thickness parameter 
= maximum camber parameter 
= camber location parameter 
= pressure coefficient 
= drag coefficient 
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Po) 

^(>7, PS) 
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Ti 


= lift coefficient 
= skin friction coefficient 
— force coefficients in x and y directions 
= camber line coordinates 
= section thickness distribution 
= surface pressure 

= velocity components in x , y and z directions 
= density 

= energy per unit volume 
= free-stream Mach number 
= free-stream density 
= free-stream velocity 
= free-stream Reynolds number 
= computational coordinates 
= horizontal blending function 
= normal blending function 
= surface slope 
= angle of attack 
= surface skin friction 


1. INTRODUCTION 


Over the past several years, Computational Fluid Dynamics (CFD) has 
rapidly evolved. This has been the result of the immense advances in computational 
algorithms [1] and the development of supercomputers. These innovations have had a 
major impact on obtaining numerical flow simulations about complex geometries. On 
current supercomputers, viscous-compressible flow simulations about aircraft config- 
uration can require several hours per steady-state solution [2]. Such large amounts 
of computational time are acceptable for proof-of-concept studies and selective anal- 
ysis, but they are not acceptable for optimization and design. With advent of the 
next generation of parallel supercomputers [3], airplane design and optimization using 
nonlinear CFD (Euler and Navier-Stokes equations) should become routine. For an 
individual component, such as a wing, it is now reasonable to consider design and 
optimization in conjunction with nonlinear CFD [4]. 

An essential element in design and optimization is acquiring the sensitivity 
of functions of CFD solutions with respect to control parameters. For aerodynamic 
surfaces, the control parameters specify the shapes of the surfaces. This affects the 
surface grid and the field grid which, in turn, affects the flow- field solution. There are 
two basic components in obtaining aerodynamic sensitivity. They are: (1) obtaining 
the sensitivity of the governing equations with respect to the state variables; and (2) 
obtaining the sensitivity of the grid with respect to the defining parameters. The 
sensitivity of the state variables with respect to the defining parameters are described 
by a linear-algebraic relation [5]. This study concentrates on the grid sensitivity and 
parameterization of aerodynamic surfaces. 
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The simplest method for calculating grid sensitivity is based upon finite- 
difference approximation. In this approach, a design parameter is perturbed from 
the nominal value, a new grid is obtained, and the difference between the new and 
the old grid is used to obtain the grid sensitivity derivatives. This direct, or brute 
force technique, has the disadvantages of being computationally intensive. It is the 
objective of the analytical approach to avoid the time consuming and costly numerical 
differentiation. In addition, the analytic derivatives are exact instead of approximate. 

Taylor et al. [6] set the stage for the development of the technique used in 
this study. For a steady-state solution of the Euler or Navier-Stokes equations, the 
sensitivity of a function of the solution with respect to the control parameters is to 
be found. The problem includes the determination of sensitivity of the grid with 
respect to control parameters. For grid generation, algebraic transfinite interpolation 
[7] is ideally suited for the study of parameterization. Parameters are sub-grouped 
according to their purposes (grid spacing control and surface shape control). The 
objective is to cast the surface parameterization in terms of design parameters rather 
than geometric variables. Geometric parameterization has only local control and 
consequently requires large number of parameters to define a grid. A specialized ap- 
plication of transfinite interpolation [8] is cast in terms of design parameters for a 
class of wing sections and wings. The grid sensitivity of this system is discussed. 


2. ALGEBRAIC GRID GF.NF.R ATTOIST 


2.1 Basic Formulation 

Structured algebraic grid generation techniques can be thought of as trans- 
formations from a rectangular computational domain to an arbitrarily-shaped physi- 
cal domain [7]. The transformations are governed by the vector of control parameters 
P. That is, 


X(£,»y,C,P) = { *(£,»?, C,P) y(Z,y,(, P) *(£,»/, C>P) } T (2.1) 

where 

0 < £ < 1, 0 < T) < 1, and 0 < £ < 1. 

A discrete subset of the vector- valued function X(£, r t jj, ( k , P) = X { x y z }J. k = X* 
is a structured grid for & = where i = 1,2, 3 • • • , L, j = 

1,2,3,- - - , Af and k = 1, 2, 3, • • • , N. 

The dominant algebraic approach for grid generation is the Transfinite In- 
terpolation scheme. The general methodology was first described by Gordon [9] , and 
there have been numerous variations applied to particular problems. The methodol- 
ogy can be presented as recursive formulas composed of univariate interpolations [10] 
or as the Boolean sum of univariate interpolations [7]. Here, we follow the Boolean 
sum representation but; for brevity, we restrict the process to two dimensions and 
omit some of the details that can be found in Ref. [7], Also, to be consistent with 
the example considered below, the parameterization is restricted to functions and first 
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derivatives at the boundaries (£ = 0, 1) and (rj = 0,1) and control in the interpolation 
functions. The transformation is 


X(£, 77 , P) = uev = u + v-uv 


(2.2) 


where 


and 


The term UV (not expanded here) is the tenser product of the two univariate interpo- 
lations. The boundary curves and their derivatives ( - and J = 

1,2 m,n = 0,1) are blended into the interior of the physical domain by the in- 
terpolation functions Q/(<f,Pg) and f3j(r/, P^). The boundary grid, the derivatives 
at the boundary grid and the spacing between points are governed by the param- 
eters |Pj P((j . The interpolation functions are controlled with the parameters 
|Po P^} • The entire set of control parameters can be thought of as a vector 

{pj p; pS p;} r = p. 


o = EE»;fePS) a ' X| ^’ P!) 


(2.3) 


/= 1 n =0 


2 1 


v = ee^(’?.p;) 3 " ,X( ^T p<j) . 


(2.4) 


J — 1 m=0 


2.2 Grid Algorithm 

An interactive univariate version of Eq.(2.2) using only the normal interpolant V is 
developed. This, known as Hermite Cubic Interpolation , matches both the function 
and its first derivative at two boundaries. An analytical approximation of the physical 
coordinates for a class of wing sections can be expressed as 


x = x.(r, P «)/?”((, PJ) + fl(r) a * l( ^ P ' ) /?i(l, Po) 


+ x ! (3,P < ,)$((,P;)+S( S ) 


0 * 2 ( 3 , 11 ) nl 




dt 


(2-5) 


0 


y = St (r,p<)^(f,pj) + fi(r) 3!, ' ( ^ p ' ) ,a;(i,p;) 

+ !/ 2 ( 5 ,P £ 2 ).3;(i.P 0 ') + 5( s ) ay3( ^ P ^ gj(t,p;) (2.6) 

where 

#(f,P3) = 2i 3 -3f 2 + l, 

0}(t, P2) = t 3 -2t 2 + t, 

$(t, P v o) = -2t 3 + 3t 2 , 
fl(t,P2) = t 3 -t 

and 

0 < t < 1. 

Five functions r = /i(0, s = f 2 { 0, A(r) = A',/ 3 (0, 5(a) = A' 2 / 4 (£) and 
* — /s(»7) an< ^ their implied defining parameters control the grid spacing on the bound- 
aries and the interior grid. Functions r and s define the grid spacing for lower and 
upper boundaries respectively, while R(r) and S(s) specify the magnitude of orthog- 
onality along those boundaries. The parameter t defines the grid distribution for the 
connecting curves between the two boundary. The quantities K x and I< 2 are param- 
eters that scale the magnitude of the orthogonality at the boundaries. Increasing K x 
and A 2 extends the orthogonality of the grid into the interior domain. Excessively 
large values of K x and K 2 can also cause the grid lines to intersect themselves, which 
is not a desirable phenomena. 

A discrete uniform distribution of the computational coordinate, can be 
mapped into an arbitrary distribution of the physical coordinate, using the specified 
control function. The essence of mapping, is that the abscissa corresponds to the per- 
centage of grid points and the ordinate corresponds to a particular control function 
which, in turn, relates to the geometric definition of the physical domain. The control 
function, can be either a specified analytical function, or for more general purposes, 


6 


a cubic- spline function. For example, the function 

e Kv - 1 

= ~pn j (' 2 ' 7 ) 

wouhl concentrate grid points close to the bottom or the top boundary depending on 
the magnitude and the sign of the constant K . Figure 2.1 shows a unit square used as 
a control domain for grid spacing. Figure 2.2 presents the parametric representation 
of the boundaries and the cubic connecting function of Eqs.(2.5) and (2.6). Appendix 
A provides a complete listing of the FORTRAN source module for this type of grid 
generation algorithm. 



Fig. 2.1 Unit square control domain Fig. 2.2 Cubic connecting function 



3. PROBLEM STATEMENT AND METHOD 

OF SOLUTION 


3.1 Theoretical Formulation 

Let Q be a solution of the Euler or Navier-Stokes equations on the domain X and 
Q* be a discrete solution on the grid X’ where 


Q* = { P P u pv pw pe X* = { x y z (3.1) 

and 

^ = R(Q-(X-),X-) = 0. (3.2) 

Here, R(Q"(X’),x*) is the residual of a steady-state solution as t — > oo. Let P be a 
vector of parameters that controls the grid X* such that X* = X(&, ijj, £fc, P) where 
£,rj and ( are computational coordinates [7]. The numerical sensitivity of a function 
F(Q(X)) with respect to the control parameters is 


The fundamental sensitivity equation containing j and described by Taylor 

et al. [6] is 


pR(Q*(X*),X*)l 

f aQ-(x-) 

l. 

pR(Q*(X*), X*)l 

pxq 

aQ-(x-) j 

1 dP 

r 

t)X* J 

\dpj 


dP } ~ °- 


(3.4) 


It is important to notice that Eq.(3.4) is a set of linear, algebraic equations 

and f d R (Q*(X )-X ) | are we u understood [6]. The 


, and the matrices g^*l(X*i') X ^ 
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quantitiy { ^ } is the solution to Eq.(3.4) given the sensitivity of the grid with 

respect to the parameters. Therefore, the grid sensitivity problem is described bv 

= GridSensitivity 

3.2 Surface Parameterization 

In Eq.(2.2), the parameterization is on the boundaries and in the interpo- 
lation functions. The most general parameterization of the boundaries would be to 
specify every grid point Xj j (i.e., each boundary grid point is a parameter). This 
conceivably could be desirable for the boundaries corresponding to an airplane surface 
to allow a design procedure to have the greatest possible flexibility. This, however, 
is impractical from a computational point of view. A compromise is to specify the 
knots of a spline function or the control polygon for a Bezier function. Even with this 
compromise, it could require hundreds of parameters for a wing. Here we propose a 
quasi-analytical parameterization in terms of design variables. For instance, a class of 
wing sections is specified by two camber-line parameters and a thickness distribution 
parameter; a wing is specified by several wing sections; and the wing surface is inter- 
polated from the sections. In this manner, an airplane component can be specified 
by tens of parameters instead of hundreds or thousands of parameters. The disad- 
vantage is that a great deal of generality is not available, but the generality is a moot 
point if computational capacity cannot accommodate it. For design and optimization 
with CFD, at the present time, it is advocated here to use a small number of design 
parameters for boundary definition. 




3.3 Grid Sensitivity 

As it is stated in the introduction, the simplest way to obtain grid sensitivity 
with respect to the parameterization is to vary the parameters and finite difference the 
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results. This is computationally expensive compared to analytically differentiating 
Eq.(2.2). Therefore, we propose the latter. Grid sensitivity can be used for grid 
adaptation, or it can be used for boundary design. For adaptation, the grid sensitivity 
with respect to those parameters that control the grid spacing and the shape of the 
field grid away from fixed boundaries are desired. The sensitivity is used to improve 
some grid-quality function of the solution. For design and optimization the sensitivity 
of the grid with respect to the parameters that define the shape of a boundary is 
desired. The sensitivity is used to improve a design function of the solution. 



4. WING-SECTION EXAMPLE 


4.1 Wing-Section Parameterization 


Much research has been devoted to the development and representation of 
wing sections. The NACA four-digit wing sections are examined for grid-generation 
parameterization. Families of wing sections are described by combining a mean line 
and a thickness distribution. The resultant expressions possess the necessary features 
that suit the problem, mainly the concise description of a wing section in terms of 
several design parameters. Reference 11 provides the general equations which define 
a mean line and a thickness distribution about the mean line. The design param- 
eters are: T = the maximum thickness, M = the maximum ordinate of the mean 
line or camber, and C = chordwise position of maximum ordinate. The numbering 
system for NACA four-digit wing-section is based on the geometry of the section. 
The first and second integers represent M and C respectively, while the third and 
fourth integers represent T. Symmetrical sections are designated by zeros for the 
first two integers, as in the case of NACA 0012 wing-section. Figure 4.1 provides a 
schematic of the section definition. The ^-coordinate is first mapped into the chord 
line x = x(r) = x(/ i(£)) forward and then reversed to cover both the top and bottom 
of the section. The mean line equation is 


M 

y c {x) = -^{2Cx - X 2 ) , x < C 
(1 - 2C + 2Cx - x 2 ) 

yc(i) = M[ - 2 , *>C. 


(4.1) 


(4.2) 
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The section thickness is given by 

T 

Vt(x) ) = —(0.2969x3 - 0.126x - 0.3516x 2 + 0.2843x 3 - 0.1015x 4 ). (4.3) 

The section coordinates are 


x,(r, Pj) = i Vi{r,P\) = y c (x) ± y T (i)- (4.4) 

Figures 4.2 and 4.3 show sample grids for NACA 0012 and NACA 8512 
section using this procedure. The orthogonality at the far-field boundary is ignored. 
For solid boundary, the orthognality is obtained using the components of unit normal 
vector at the surface 

dx,(r,P<)/ ' 

(4.5) 

Figure 4.4 shows a wing-surface grid derived from three differently-specified 
sections in the spanwise direction. The surface grid results from the distribution 
function/i(£ ), and interpolation of the design parameters for the three wing sections 
in the spanwise direction. In addition to the three design parameters for each wing 
section, it is necessary to specify their relative chord lengths and positions. The 
additional design parameters can be: leading edge sweep, trailing edge sweep, dihedral 
angle, reference chord length, and section-span locations. 


ari(r,Pj) 

dt 


= 4= stnd 


cMr, Pj| 

dt 


— ±COS0 


0 = tan' 


Leading edge 



Fig. 4.1 Schematic of wing-section 




Fig. 4.2 Example grid for NACA 0012 wing-section 
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Fig. 4.3 Example grid for NACA 8512 wing-section 
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Fig. 4.4 Example wing surface grid (top view) 

Appendix B provides the source module for generating the surface of NACA four-digit 
wing-sections. 

4.2 Grid Sensitivity With Respect To Control 

Parameters 

There are two types of control parameters involved in this analysis. First, there 
are the design parameters (T, M, and C) which specify the shape of the piimai} 
boundary and secondly, there are the parameters that define the other boundaries 
and the the spacing between grid points. Here we express, in part, the sensitivity of 
the grid with respect to the design parameter vector Xp, and with k a stretching 
parameter in the interpolation functions related by / 5 (r/,k)). The grid sensitivity 
with respect to design parameters at the outer boundary has been ignored. Also, due 
to zero orthogonality at the outer boundary, a direct differentiation of Eqs.(2.5) and 
(2.6) with respect to Xp yields to 

dx dx dxi(r, Pi) dx dx\(r, P}) 

dX n = dTt dX D + dx\ dX D 


(4.6) 


14 


_dy_ = dy dy^Pj) d^dy^ P|) 

dx D dy[ dX D + d Vl dx D 

where 


X D = (T, M, C) 


dx dy n , 

a^ = ^ = M) 


dx dy . 

l)x{ = ~dy[ = R M*l). 


;i-8) 

11-9) 


The prime indicates differentiation with respect to t and can be substituded from 
Eq.(4.5). Since x^r, Pj) is independent of design parameters Xp, then 

(r, PS) 


dX D 


= 0 . 0 . 


( 1 . 10 ) 


The x coordinate sensitivity, Eq.(4.6), can now be reduced to 

0 = tan 

- ” CJA D 

Using the relation 




-1 l dyi(r,pj) ' 


d 


tan l u = 


dX D 

the x coordinate sensitivity becomes 

dx 


1 


du 


1 + u 2 dX D 


dX 


D 




d d yi (r,P[) 


' eni^ ydXydx^Pj) 


1.11 


( 1 - 12 ) 


(1.13) 


1 + 


The term can be evaluated by direct differentiation of Eq.(4.4). The y 

coordinate sensitivity with respect to design papameters can be obtained using similar 
procedure. Equation (4.7) can be modified to 

d d yi (r,P\ 


fijr- = P o) ~ ] P5)«'ng- 


1 


1 + 


" a Vl (r,P^) \ 2 dX D dxi{r, Pf ) 
ax,(r,p< ) J 

(1.14) 

All terms at the right hand side of Eqs.(4.13 ) and (4.14 ) can be evaluated explicitly 
due to analytical parameterization of the surface for this particular example. The 
grid sensitivity with respect to the stretching parameter k are 

a', rP < 1 «’?(‘.P!)« &:,(>•, p«)3/j}(<,P3)a< 

ak “ Il(r ' Pl, ~ ak + R{r) — Wt — ak 
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+ x 2 (r,P^ 


dd%(t, PS) dt 


dt 


dk 


r + S(s) 


dx 2 (s,Pl)d0](t,P' o ) dt 


dr) 


dt dk 


( 4 . 15 ) 


where 


aftt, p;, 

dt 

Mik p o) 




= 6r - 6< 

= — 6t 2 6 f 


d$Ut,P v 0 ) „ 

- u - 0 = 3< 2 —44 + 1 


dt 

mu pg) 

dt 


= 3 r - l. 


An example of the stretching function is 


e k,) - 1 

t = r 

e k - 1 

d* _ (e k - l)i/e £,> - (e kr? - l)e k 
dk ~ (e k - l)* 


(4.16) 

(4.17) 


Similar developments can be extended to other grid control parameters such 
as the distribution of grid point around the wing section and magnitude of orthogo- 
nality at the boundaries. Appendix C provides the source module for grid sensitivity 
of NACA four-digit wing-sections with respect to design parameters. 


4.3 Flow Sensitivity With Respect To Control 

Parameters 

The flow sensitivity coefficient { ^ } can now be evaluated using the 

fundamental sensitivity equation. Eq.(3.4). The Jacobian matrices 

j r<9R(Q*fx*) x*n 

an d ax* — 1 are °btained using an implicit time integration of the 2-D thin- 

layer Navier-Stokes equations [10]. These equations are solved in their conservation 
form using an upwind cell-centered finite- volume formulation. A detailed description 
of the procedure is found in Refs. [12-17] and is not repeated here. A third-order 
accurate upwind biased inviscid flux balance is used in both streamwise and normal 
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directions. The finite- volume equivalent of second-order accurate central differences is 
used for viscous terms. For a typical design analysis of an airfoil, the flow sensitivity 
coefficient j provides far more information than needed. In most cases, the 

sensitivity of aerodynamic forces on the surface, such as lift and drag coefficients, 
are sought. For such analysis, only a small subset of the flow sensitivity coefficient 
(j e sur f ace properties) is needed since the lift and drag coefficients can 
be expressed as 

Cl = Cycosa — Cxsina (4.18) 


Cd = Cysina + Cxcosa 


(4.19) 


where a is the flow angle of atack. The quantities Cx and Cy are the total force 
coefficents along x and y directions respectively and can be expressed as 


NE 

Cx = ^2C Pt {y i+ i - y t ) + C fi (xi+i -*,) 

t=l 

(4.20) 

.ME 

Cy = £C Pl (x, +1 -x,-) + Cj l {y i+1 -y t ) 

(4.21) 


1=1 


where C pt and C j t are pressure and skin friction coefficients respectively 


c„ 


Pi 

\pooUlo 


Cu = 


\pooUl, 


(4.22) 


and N E represents total number of bondary cells along airfoil surface. The terms P, 
and T, are pressure and shear stress associated with boundary cell i and the quantity 
\ Poo Ul is known as dynamic pressure of the free stream . Finally, the drag and lift 
sensitivity coefficients with respect to Xd are obtained by differentiating Eqs.(4.18) 
and (4.19) as 


dC L 

dC Y 

dCx 

(4.23) 


= dx u COSQ 

— sina 

dXn 

dC L 

3Cy 

dC x 

(4.24) 

<9X d 

= , w cosa - 
aX D 

sma. 

oX D 



5. RESULTS AND DISCUSSION 


5.1 NACA 0012 Airfoil Test Case 

5.1.1 Grid Sensitivity 

The first test case considered is the NACA 0012 symmetrical airfoil. The 
previously obtained grid, as shown in Fig. 4.2, is considered for grid sensitivity anal- 
ysis. The grid sensitivity with respect to the vector of design parameters X D . are 
obtained using Eqs. (4.13) and (4.14). The maximum thickness T is the only design 
parameter for this case. 

Figure 5.1 shows the contour levels of the y-coordinate sensitivity with re- 
spect to the thickness parameter, T. The highest contour levels are, understandably, 
located in the vicinity of the chordwise location for the maximum thickness of the 
wing section. For a NACA four-digit wing section, this is positioned about 0.3 of 
the chord length from the leading edge [11]. The positive and negative contour levels 
corresponding to the upper and lower surfaces are the direct consequence of Eq.(4.4) 
and the second term on the right hand side of Eq.(4.7). The sensitivity levels decrease 
when approaching the far-field boundary due to diminishing effects of the interpola- 
tion function /?°(t, P^). The first term on the right hand side of Eq.(4.7) is responsible 
for the sensitivity effects due to orthogonality on the surface, and it is directly pro- 
portional to the magnitude of the orthogonality vector K \ . The wake region is not 
sufficiently affected by any of the design parameters, and no major sensitivity gradi- 
ent should be expected there. 
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Figure 5.2 demonstrates the contour levels of the x-coordinate sensitivity 
with respect to thickness parameter, T. An interesting observation can be made here 
regarding the contour levels adjacent to the surface. Unlike the y-coordinate sensitiv- 
ity, the x-coordinate sensitivity has its lowest value on the surface. This can be traced 
back to Eq.(4.4), which indicates that x-coordinates on the surface are basically inde- 
pendent of the design parameters. The only remaining factor is the second term on 
the right hand side of Eq.(4.6) which is the effect of the surface orthogonality vector. 
There are some negative pockets of contour levels on the forward section and some 
corresponding positive pockets on the rear section. The dividing line between these 
pockets is located near the location of the maximum thickness (i.e., 0.3 of chord from 
leading edge). A simple conclusion from Fig. 5.2 is that by increasing the thickness 
parameter, T points on the forward section will move to the left, while at the same 
time, points at rear section will move to the right. 

For comparison purposes, the grid sensitivity for this case is obtained using 
the finite difference approach. The design parameter(i.e. T for this case) are per- 
turbed, one at a time, and a new grid is obtained using Eqs.( 2.5) and (2.6). The 
sensitivity is then computed using a central difference approximation and the results 
are presented in Figs. 5.3 and 5.4 . A side by side comparison of both results indicates 
good agreement between the two approaches. 

5.1.2 Flow Sensitivity 

The second phase of the problem is obtaining the flow sensitivity coefficients 
using the previously obtained grid sensitivity coefficients. In order to achieve this, 
according to Eq. (3.4), a converged flow field solution about a fixed design point 
should be obtained. The computation is performed on a C-type grid composed of 141 
points in the streamwise direction and 31 points in the normal direction. It is appar- 
ent that such a coarse grid is inadequate for capturing the full physics of the viscous 
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flow over an airfoil. Therefore, it should be understood that the main objective here 
is not to produce a highly accurate flow field solution rather than to demonstrate the 
feasibility of the approach. 

A free stream Mach number of M ^ = 0.8, Reynolds number Re^, = 10 6 . 
and angle of attack a — 0° is used. Figures 5.5 and 5.6 demonstrate the pressure and 
Mach number contours of the converged solution. Figure 5.7 shows the surface pres- 
sure coefficient C p , where the lift and drag coefficient for this particular example are 
Cl = 1.53x10 , Cd = 4.82xl0~ 2 . The sensitivities of the aerodynamic forces, such 

as drag and lift coefficients with respect to thickness parameter T, are obtained uti- 
lizing Eqs. (4. 18-4.24). The corresponding results are presented in Table 5.1. Again, 
for comparison purposes, a finite difference approximation has been implemented to 
validate the results. A nominal perturbation of 10 -3 for design parameter T has been 
chosen and the corresponding results are included in Table 5.1. The good agreement 
between the two sets of numbers verifies the accuracy of the approach. 

Another important goal of using sensitivity analysis, apart from optimiza- 
tion, is the approximation analysis. An approximate version of Eq.(3.4) can be used 
to predict the steady-state solution changes which occur in response to geometric 
shape changes. Such a method is valid as long as the changes in geometric shape 
(i.e., design parameter) are small. Figure 5.8 shows the non-linear relation between 
drag coefficient and thickness parameter T verifying the above argument. 


5.2 NACA 8512 Airfoil Test Case 

5.2.1 Grid Sensitivity 

The second test case considered is the NACA 8512 cambered airfoil. Again, 
the previously obtained grid, as shown in Fig. 4.3, is considered for grid sensitivity 
analysis. Figures 5.9 and 5.10 show the y and x-coordinate sensitivity with respect 
to parameter T respectively. Their characteristics are similar to the previous sym- 
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metrical airfoil case; hence, detailed description of their behavior is omitted here. 

Figure 5.11 represents the y-coordinate sensitivity with respect to camber, 
M. It appears that the highest sensitivity contour levels are located at the chordwise 
location of camber, C (i.e., 0.5 of chord length). The contour levels decrease toward 
the far-held boundary, again as a consequence of interpolation function. However, 
unlike Fig. 5.9 here they possess positive values on both upper and lower surfaces. 
Consequently, an increase in camber, M, shifts all the points upward. Again, mini- 
mum activity can be detected in the wake region. Figure 5.12 shows the x-coordinate 
sensitivity contours with respect to camber, M. Here, as in Fig. 5.10, the sensitivities 
are minimum on the surface of the wing-section. There is a small gradient on the 
forward section, but by far, the strongest gradient is in the rearward section due to 
orthogonality effects. 

Figure 5.13 illustrates the y-coordinate sensitivity with respect to camber 
location, C. A dividing line between positive and negative contour levels appears 
near the chordwise position of the camber. Like previous cases, there is no significant 
activity in the w'ake region. The result indicates that a positive change of C will cause 
the movement of points downward on the forward section, while at the same time, the 
points on the rear section will respond by moving upward. Figure 5.14 illustrates the 
x-coordinate sensitivity with respect to camber location C. The two major features 
are attributed to chordwise location of the camber and the orthogonality effects on 
the tail section. It is interesting to notice that the sensitivity level for camber location 
is considerably less than the other two design parameters. 

Similar developments can be extended to other grid control parameters such 
as the distribution of grid point around the wing section and magnitude of orthogo- 
nality at the boundaries. For example, the grid sensitivity with respect to stretching 
parameter k, using Eqs.(4.15-4.17), is obtained and the results are presented on Figs. 
5.15 and 5.16. 



5.2.2 Flow Sensitivity 

Using free stream conditions of M ^ = 0.8 , Re ^ = 10 6 , and a = 0°, a 
converged flow field solution is obtained. As in previous case, a C-type grid of 141x31 
is used. Figures 5.17 and 5.18 illustrate the pressure and Mach number contours. 
Figure 5.19 shows the surface pressure coefficient C p , where lift and drag coefficients 
are Cl — 0.106, and Co = 0.0738. The aerodynamic sensitivity coefficients with 
respect to vector of design parameters Xd, are obtained and presented in Table 5.2. 
A comparison with finite difference validates the accuracy of the approach. 
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Fig. 5.1 Y-coordinate sensitivity with respect to thickness T (Analytical) 
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Fig. 5.2 X-coordinate sensitivity with respect to thickness T (Analytical) 




21 




1.06457 
1 .02346 
0.982351 
0.941241 
0,900131 
0.859021 
0.817911 
0.776801 
0.735691 
0.694582 
0 653472 
0.612362 
0.571252 
0.530142 
0.489032 


Flow Condition: 
Angle of Attack = 0.0 
Mach Number = 0.8 
Re = 1000,000 


Fig. 5.5 Pressure contours for NACA 0012 wing-section 
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Fig. 5.6 Mach number contours for NACA 0012 wing-section 




ig. . r ,.7 Surface pressure coefficient for NACA 0012 wing-section 
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Fig. 5.11 Y-coordinate sensitivity with respect to camber M (Analytical) 
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Fig. 5.12 X-coordinate sensitivity with respect to camber M (Analytical) 
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Fig. 5.13 Y-coordinate sensitivity with respect to camber location 
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Fig. 5.14 X-coordinate sensitivity with respect to camber location C (Analytical) 





Fig. 5.17 Pressure contours for NACA 8512 wing-section 



Fig. 5.18 Mach number contours for NACA 8512 wing-section 
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Fig. 5.19 Surface pressure coefficient for NACA 8512 wing-section 
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6. CONCLUSION 


The objective of this study has been to demonstrate an approach for ob- 
taining grid sensitivity which can be used in aerodynamic design and optimization. 
It is shown that grid sensitivity is an essential ingredient in the calculation of aero- 
dynamic sensitivity. The main supposition is that a grid is defined algebraically in 
terms of parameters and computational coordinates. Therefore, coordinates of the 
grid and derivatives of the coordinates with respect to the parameters (grid sensi- 
tivity) are computed directly as functions of the parameters and uniformly-discrete 
values of the computational coordinates. A subset of the parameters defines the shape 
of the grid boundaries which corresponds to the aerodynamic surfaces of interest. It 
is recommended that the aerodynamic surfaces be parameterized in terms of design 
parameters which have global control. As compared to a geometric parameteriza- 
tion, this drastically reduces the number of parameters. However, it limits the design 
flexibility. In addition to the aerodynamic surface parameters, the sensitivity with 
respect to parameters that define other boundaries, such as the far-field boundary or 
the spacing of grid points, is available for analysis or grid adaptation. 

The algebraic grid-generation scheme and NACA wing sections presented 
here are intended to demonstrate the elements involved in obtaining grid sensitivity 
from an algebraic grid generation process. It is evident that each grid generation 
formulation would require considerable analytical differentiation. This implies that 
a symbolic manipulator, which directly produces computer code for derivative eval- 
uation, should be considered. Also, there are trade-offs between analytical differen- 
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tiation and finite-difference differentiation. It may be feasible to obtain some of the 
derivatives by finite differences. 

It is implied that airplane surfaces should be parameterized in terms of de- 
sign variables. This is not simple or feasible for geometrically-complex airplanes in 
advanced stages of design. However, design parameterization is feasible during con- 
ceptional and preliminary design. The parameterization, which is the development 
of analytical formulas for part or all of a surface is critical for satisfactory results. 
It is always possible to create a geometric parameterization of a surface (collection 
of points or derivatives that define surface patches), but geometric parameterizations 
have very local sensitivities and a complete aerodynamic surface can require a large 
number of parameters for its definition. 

As a compromise between totally analytic parameterization of surfaces and 
geometric parameterization, a hybrid approach is advocated. In a hybrid approach, 
certain sections or skeletal parts of a surface are specified analytically and interpo- 
lation formulas are used for the remainder of the surface. This is employed for the 
wing example described herein. 
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APPENDIX A 

FORTRAN LISTING FOR GRID GENERATION 
ALGORITHM : HERMITE 
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SUBROUTINE HERMITE1 (XS, YS, NI, NJ, ALFA, XGRD, YGRD, PROBID, DOBATCH) 
PARAMETER (NGRID=3) 

$ INCLUDE tbggl.inc 
$ INCLUDE tbggS.inc 

DIMENSION XS (NGDIM, *) , YS (NGDIM, *) , T (NGDIM) , XGRD (NGDIM, NGDIM) 

1 , YGRD (NGDIM, NGDIM) , ALFA (NGDIM) , SI (NGDIM, NGDIM) 

2 , XLEFT (NGDIM) , YLEFT (NGDIM) , XRIGHT (NGDIM) , YRIGHT (NGDIM) 

3 , RI (NGDIM) , RO (NGDIM) , R (NGDIM) , S (NGDIM) , SK (NGDIM) 

4 , X (NGDIM) , Y (NGDIM) , XVIEW (NGDIM, NGDIM) , YVIEW (NGDIM, NGDIM) 
COMMON/PQ/PS (MGDIM) , QS (MGDIM) , t S (MGDIM, MGDIM) 
COMMON/DXYDETA/DXBDETA (NGDIM) , DYBDETA (NGDIM) 

COMMON/ JAYS / J1 , J2 , J3 

common/design/ cm, p, th, rr, wlength, chord 
CHARACTER PROBID* (*) , PR (4) *80 
SAVE RO, XLEFT, YLEFT, XRIGHT, YRIGHT 
LOGICAL REDRAW, DONE, RED I ST, DOBATCH, FIRST 

C 

c 

c 

c 

C Does basic grid calculations, draws grid, and provides for user 
C modifications in interactive loop. 

C 

C 

C 


c 

REDIST 

= 

.FALSE 


REDIST 

= 

.TRUE. 

c 

DONE 

35 

.FALSE 


DONE 

35 

.TRUE. 


c 

c 

c 

c 


c 


299 


Magnitude Of Orthogonality Vector 


WRITE (*, *) 

WRITE (*, *) ' Magnitude of Normal Derivatives ? ' 
READ ( * , *) SKB 

DO 299 1=1, NI 
SK (I) = SKB 
CONTINUE 
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99 


C 

C 

C 

C 

c 

c 


c 

c 


CONTINUE 


Distribution For Stretching Variable T 
Uses Linear Distribution For Initial Trial and Arc-length 
Disribution For Final Trial 


IF (RED 1ST .AND . ( .NOT. DONE) ) THEN 

DONE = .TRUE. 

IF (DOBATCH) THEN 

CALL BATCH (NJ, RO, 5) 
ELSE 

CALL DISTXI (NJ, RO, 5) 
END IF 


C 


DO 100 I = 1 , NI 


o o 
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C 


199 

C 

C 

c 


c 

101 

c 

c 

c 


c 

102 

c 

100 

c 

c 


DO 199 J = 1 , NGDIM 
X(J) = XGRD (I, J) 

Y(J) = YGRD (I, J) 

CONTINUE 

CALL ARC (X, Y, R, NGDIM, RMAX) 

DO 101 J = 1 , NGDIM 

RI ( J) ■ R ( J) / RMAX 

IF (RI (NGDIM) .NE.1.0)RI (NGDIM) =1.0 

CONTINUE 

CALL INTERPO(RI, T, NGDIM, RO,S,NJ) 

DO 102 J = 1 , NJ 

S1(I,J) = S(J) 
tS(I,J) = S1(I,J) 

CONTINUE 


CONTINUE 

ELSEIF ( ( . NOT . REDIST) .AND. ( .NOT. DONE) ) THEN 
REDIST = .TRUE. 


C 

C 


C 

201 

C 

C . . 
C 

C 


C 


107 


C 


DO 201 1=1, NI 

DO 201 J = 1 , NGDIM 

T ( J) = FLOAT (J-l) /FLOAT (NGDIM-1) 
SI (I , J) = T ( J) 
tS(I,J) = S1(I,J) 


CONTINUE 

Base Line Distribution 

ELSEIF (REDIST .AND .DONE) THEN 

IF (DOBATCH) THEN 

CALL BATCH (NJ, RO, 5) 

ELSE 

CALL DISTXI (NJ, RO, 5) 

END IF 

DO 107 1=1, NI 
DO 107 J = 1 , N J 
S1(I,J) = RO ( J) 
tS(I,J) = S1(I,J) 

CONTINUE 

END IF 

IF (MGDIM.NE. NGDIM) THEN 

WRITE ( * , *) 9 »» tbggl.inc 6 tbggS . inc Should Have the 
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C 

C 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 


c 

c 


c 


c 

c 

c 

c 

c 


1 Same Dimension' 

STOP 
END IF 


Hermite Interpolation Function (Eq.3 of manual) 


SKT =0.0 

DO 240 1=1, NI 

XB = XS ( 1 , 1) 

YB = YS (1,1) 

XT = XS ( I , 2 ) 

YT = YS ( I, 2 ) 

Orthogonality For Foil 

IF ( X . GE . 1 . AND . X . LE . J2 ) THEN 

DXDETA = - SIN (ALFA ( I) ) 
DYDETA = + COS (ALFA (I)) 

ELSEIF ( I . GT . J2 . AND . I . LE . NI ) THEN 

DXDETA = + SIN (ALFA (I)) 
DYDETA = - COS (ALFA (I)) 


END IF 


DXBDETA ( I ) =DXDETA 
DYBDETA ( I ) =DYDETA 


IF (I .EQ . J2-4) SK (I) 

= 

SKB 

* 

0.85 

IF (I . EQ . J2-3) SK (I) 

= 

SKB 

* 

0.75 

IF (I .EQ. J2-2) SK(I) 

S 

SKB 

* 

0.65 

IF (I . EQ . J2-1) SK (I) 

= 

SKB 

* 

0.55 

IF (I . EQ . J2) SK (I) 

= 

SKB 

* 

0.45 

IF (I . EQ . J2+1) SK (I) 

= 

SKB 

* 

0.55 

IF ( I . EQ . J2+2 ) SK ( I ) 

= 

SKB 

* 

0.65 

IF (I . EQ . J2+3) SK (I) 

= 

SKB 

★ 

0.75 

IF (I . EQ . J2+4) SK (I) 

= 

SKB 

★ 

0.85 


PS(I) = SK (I) 

QS (I) = SKT 

IF ( . NOT . DONE) THEN 

N=NGDIM 

ELSE 

N=NJ 

END IF 


C 


DO 240 J = 1 , N 



o o 
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C 

C 


C 

240 

C 


CALL BLENDF (SI (I, J) ,F1,F2,F3,F4) 

XGRD (I, J) = XB*Fl + XT*F2 + PS (I) * DXDETA 
YGRD ( I , J) = YB*F1 + YT*F2 + PS (I) * DYDETA 


CONTINUE 

IF ( .NOT. DONE) GOTO 99 


RETURN 

END 



APPENDIX B 

FORTRAN LISTING FOR NACA FOUR-DIGIT SURFACE 

GENERATION : NACA 
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SUBROUTINE NACA (IS, XS, YS, NLI, ALFA, DOBATCH, CHEAT) ^ 

$ INCLUDE tbggl.inc 

PARAMETER (NGRID=3) 

DIMENSION XS (NGDIM, 4) , YS (NGDIM, 4) ,NLI (4) , ALFA (NGDIM) 

1 , XOU (NGDIM) , YOU (NGDIM) , XOL (NGDIM) , YOL (NGDIM) 

2 , YC (NGDIM) , XCI (NGDIM) , YT (NGDIM) , DYTDXCI (NGDIM) 

3 , DYCDXCI (NGDIM) , ETA (NGDIM) 

DIMENSION XI (NGDIM) , YI (NGDIM) , XO (NGDIM) , YO (NGDIM) , RI (NGDIM) 

1 , RO (NGDIM) , ROBAR (NGDIM) , X (NGDIM) , XX (NGDIM) 

2 , XU (NGDIM) , YU (NGDIM) , XL (NGDIM) , YL (NGDIM) , XDIST (NGDIM) 

COMMON / A/ A1 , A2 , A3 , A4 , A5 

COMMON /DESIGN/CM, P , TH, R, WLENGTH, CHORD 

COMMON / TANGENT / ALF AU (NGDIM) , ALFAL (NGDIM) , DYUDX (NGDIM) , DYLDX (NGDIM) ,NFOIL 
COMMON/OLD/XUS (NGDIM) , YUS (NGDIM) , XLS (NGDIM) , YLS (NGDIM) 

LOGICAL RED 1ST (3) , DOBATCH , CHEAT 
DATA PI,CX/3.14159,7./ 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


»» NACA FOUR DIGIT AIRFOIL SECTION «« 

COMPUTES AN AIRFOIL SECTION ANALITICALLY USING THE 
RELATIONS GIVEN IN 'THEORY OF WING SECTIONS' 


WRITTEN BY : IDEEN S ADREHAGH IGH I 

MECHANICAL ENGINEERING & MECHANICS 
OLD DOMINION UNIVERSITY 
7-12-1991 


INPUT PARAMETERS 


CHORD =1.0 

WRITE (*,*)' r ,= == = = r- ' 

WRITE ( * , *) ' »»» AIRFOIL GEOMETRY «««' 

WRITE (*, *) ' =======--— TT — = 

WRITE (*, *) 

WRITE (*,*) 'NUMBER OF POINTS IN XI -DIRECT ION? ' 

READ (*, *)NI 
NFOIL=NI 

WRITE (*,*) 'NUMBER OF POINTS IN ETA-DIRECTION?' 

READ ( * , * ) N J 

WRITE (*,*) 'OUTER BOUNDARY LOCATION (CHORD LENGTH)?' 

READ (* , *) R 
R=CHORD*R 
YLENGTH=R 
PERCENT =CHORD / 1 0 0 . 

TENTH=CHORD / 1 0 . 

WRITE (*,*) 'MAXIMUM ORDINATE OF MEAN LINE OR CAMBER (% CHORD) ?' 
READ ( *, *) CM 
CM=CM* PERCENT 

WRITE (*,*)' CHORDWISE POSITION OF CAMBER (./CHORD)?' 

READ (*, *) P 
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C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

c 

c 

c 


c 

c 

c 

c 


5 

c 

c 

c 


c 

c 

c 


p=p * tenth 

WRITE (*, *) ' MAXIMUM THICKNESS OF AIRFOIL ( % CHORD) ?' 

READ (*, *) TH 

TH=TH*PERCENT 


AIRFOIL NOMENCLATURE: 


CM CAMBER 

p CAMBER LOCATION ALONG CHORD 

TH MAXIMUM THICKNESS 

A THICKNESS DISTRIBUTION COEFFICIENTS 

XCI , YC MEAN LINE COORDINATES 

CHORD CHORD LENGTH 

WLENGTH WAKE LENGTH 


- THICKNESS DISTRIBUTION 

A1=0. 29648 
A2=-0 . 12642 
A3=-0. 35202 
A4=0. 28388 
A5=-0 . 10192 


INITIAL DISTRIBUTION IN XI-DIRECTION FOR LOWER BOUNDARY 

DO 5 1=1, NGDIM 

C = FLOAT ( I -1) /FLOAT (NGDIM- 1) 

X(I) = (EXP (CK*C) -1 . ) / (EXP (CK) -1 . ) 

XCI (I) = X(I) * CHORD 
CONTINUE 

GET INITIAL FOIL (UPPER £ LOWER) 

CALL FOIL (XCI, XUS, YUS, ALFAU,DYUDX, NGDIM, 1, CHEAT, 999) 

CALL FOIL (XCI, XLS,YLS,ALFAL,DYLDX, NGDIM, 2, CHEAT, 999) 

OUTER BOUNDARY 


C 

C 

C 

C 

c 


c 

35 

C 


XOUTER = CHORD+R 


• • • • INITIAL DISTRIBUTION IN XI-DIRECTION FOR OUTER BOUNDARY 
DO 35 1=1, NI 

C = FLOAT ( I -1) /FLOAT (NI-1) 

XX (I) = (EXP (C*CK) -1 . ) / (EXP (CK) -1 . ) 

CONTINUE 

DO 36 1=1, NI 

XOL(I) = XX(I) * XOUTER 


44 


C 


XOU ( I ) =XOL ( I ) 
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IF(XOL(I) . LE.R) THEN 

YOU (I) =SQRT (R*R- (XOL (I) -R) * (XOL (I) -R) ) 
YOL (I) =-YOU (I) 

ELSE 

YOU ( I ) =YLENGTH 
YOL ( I ) =- YLENGTH 
END IF 


36 CONTINUE 

C 

C OUTPUT TO HERMITE 

C 

C NLI ( J) = REPRESENTS NUMBER OF POINTS IN EACH BOUNDARY 

C J=1 SOLID BOUNDARY 

C J=2 OUTER BOUNDARY 

C J=3 RIGHT BOUNDARY 

C J=4 LEFT BOUNDARY 

C 


NLI (1) = NI*2 
NLI (2) = NI*2 
NLI (3) = NJ 
NLI (4) = NJ 
C 

C BOTTOM BOUNDARY DISTRIBUTION (ARC-LENGTH) 

C 

J=1 

C 

C REDIST(J) = .TRUE. 

REDIST(J) = .FALSE. 

C 

IF (REDIST ( J) ) THEN 
C 

C INTERPOLATE UPPER PORTION 

C 

CALL ARC (XUS, YUS, RI, NGDIM, RMAX) 

IF (DOBATCH) THEN 

CALL BATCH (NI , ROBAR, 1 ) 

ELSE 

CALL DISTXI(NI, ROBAR, IS) 

END IF 

DO 100 1=1, NI 

RO (I) = ROBAR (I) * RMAX 


100 CONTINUE 

CALL INTERPO (RI, XUS, NGDIM, RO, XDIST, NI) 

C 

C GET FINAL FOIL (UPPER) 

C 

CALL FOIL (XDIST, XU, YU, ALFAU, DYUDX, NI, 1, CHEAT, 999) 
C 

C INTERPOLATE LOWER PORTION 

C 


CALL ARC (XLS,YLS,RI, NGDIM, RMAX) 
IF (DOBATCH) THEN 

CALL BATCH (NI, ROBAR, 2) 

ELSE 

CALL DISTXI (NI, ROBAR, IS) 

END IF 

DO 101 1=1, NI 

RO (I) = ROBAR ( I ) * RMAX 
101 CONTINUE 


File: a Printed Tue Jan 01 00:00:00 198 Page: 4 


C 

C 

c 

c 


CALL INTERPO (RI, XLS, NGDIM, RO, XDIST, NI) 

GET FINAL FOIL (LOWER) 

CALL FOIL (XDIST, XL, YL, ALFAL, DYLDX, NI, 2, CHEAT, 999) 
ELSE 


C 

C Base Line Distribution 

C 

IF (DOBATCH) THEN 

CALL BATCH (NI , ROBAR, 1 ) 

ELSE 

CALL DISTXI (NI, ROBAR, IS) 

END IF 

CALL FOIL (ROBAR, XU, YU, ALFAU, DYUDX, NI, 1, CHEAT, 999) 

IF (DOBATCH) THEN 

CALL BATCH (NI, ROBAR, 2) 

ELSE 

CALL DISTXI (NI, ROBAR, IS) 

END IF 

CALL FOIL (ROBAR, XL, YL, ALFAL, DYLDX, NI, 2, CHEAT, 999) 

END IF 
C 

C SHIFT X-COORDINATES TO THE RIGHT FOR C-TYPE GRID 

C 

DO 30 1=1, NI 


XU(I) = XU (I) + R 
XL (I) = XL (I) + R 
30 CONTINUE 
C 

C ASSEMBLE COUNTER-CLOCKWISE 


C 

DO 40 1=1, NLI ( J) 

IF (I . LE .NI) THEN 

XS(I, J)=XU(NI-I+1) 

YS (I, J) =YU (NI-I+1) 
ALFA ( I ) =ALFAU (NI - 1+1 ) 
ALFA(NI) =PI/2 . 

ELSE 

XS(I, J)=XL(I-NI) 

YS (I, J) =YL (I-NI) 

ALFA ( I ) =ALFAL ( I -NI ) 
END IF 


40 CONTINUE 
C 

C CHECK DOUBLE POINTS FOR BOTTOM BOUNDARY 

C 


DO 46 1=1, NLI (J) -1 

IF (XS (I , J) . EQ .XS (X+l, J) .AND . YS (I, J) .EQ . YS (1+1, J) ) THEN 
K=I+1 

XS (K, J) =XS (K+l, J) 

YS (K, J) =YS (K+l , J) 

ALFA (K) =ALFA (K+l ) 

IF (I . EQ.NLI ( J) -1)NLI ( J) =NLI ( J) -1 

END IF 


46 CONTINUE 
C 

C OUTER BOUNDARY 

C 


46 
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J=2 

C 

REDIST(J) = .TRUK. 

DO 45 1=1, NLI (J) 

IF(I.LE.NI) THEN 

XI (I) =XOU (NI-I+1) 
YI(I)=YOU(NI-I+l) 
ELSK 

XI(I)=XOL(I-NI) 

YI (I) =YOL (I-NI) 
END IF 


45 CONTINUE 
C 

C CHECK FOR DOUBLE POINTS (SYMMETRY LINE) 


C 

DO 48 1=1, NLI (J) -1 

IF (XI (I) . EQ.XI (1+1) .AND . YI (I) .EQ.YI(I+1)) THEN 
K=I+1 

XI (K) =XI (K+l) 

YI (K)=YI (K+l) 

IF (I . EQ .NLI ( J) -1) NLI ( J) =NLI ( J) -1 
END IF 

48 CONTINUE 
C 

IF (REDIST ( J) ) THEN 

CALL ARC (XI, YI, RI, NLI ( J) , RMAX) 

IF (DOBATCH) THEN 

CALL BATCH (NLI (J) , ROBAR, 6) 

ELSE 

CALL DISTXI (NLI (J), ROBAR, IS) 

END IF 

DO 103 1=1, NLI (J) 

RO (I) = ROBAR (I) * RMAX 
103 CONTINUE 

CALL INTERPO (RI, XI, NLI ( J) , RO, XO, NLI ( J) ) 

CALL INTERPO (RI, YI, NLI ( J) , RO, YO, NLI (J) ) 

END IF 
C 

DO 201 1=1, NLI (J) 

IF (REDIST ( J) ) THEN 
XS(I,J) = XO(I) 

YS(I,J) = YO (I) 

ELSEIF ( .NOT .REDIST (J) ) THEN 
XS(I,J) = XI (I) 

YS (I, J) = YI (I) 

END IF 
201 CONTINUE 


C 

C ETA DISTIBUTION 

C 

DO 51 J = 1 , NJ 

ETA(J) = FLOAT (J-l) /FLOAT (NJ-1) 
51 CONTINUE 
C 

C RIGHT BOUNDARY 

C 

J=3 


DO 49 1=1, NJ 
XS (I, J)=CHORD+R 
YS(I, J)=ETA(I)*R 
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49 CONTINUE 4 

C 

C LEFT BOUNDARY 

C 

J=4 

DO 50 1=1, NJ 
XS (I, J)=CHORD+R 
YS(I, J)=-ETA(I)*R 

50 CONTINUE 
C 

C 

RETURN 

END 

C 

Q************************************************************************ 

c 

SUBROUTINE FOIL (XCI, X, Y, ANGLE, DYDX, N, II, CHEAT, I GRID) 

$ INCLUDE tbggl.inc 

DIMENSION XCI(NGDIM) , X (NGDIM) , Y (NGDIM) , YC (NGDIM) , YT (NGDIM) ,XD (3) 

1 , DYCDXCI (NGDIM) , DYTDXCI (NGDIM) , ANGLE (NGDIM) , DYDX (NGDIM) 

COMMON /A/A1,A2, A3, A4, A5 
COMMON /DESIGN/CM, P , TH, R, WLENGTH, CHORD 
COMMON /DELTAXD/TH1, CM1, PI, TH3,CM3, P3 
LOGICAL CHEAT 
DATA PI, CK/3 . 14159, 3 . / 

C 

C 

C 

XD (1) = TH 
XD (2) = CM 
XD (3) = P 
C 

IF (CHEAT) THEN 

IF ( IGRID . EQ . 1 ) THEN 
XD (1) = TH1 
XD (2) = CM1 
XD (3) = PI 

ELSE IF ( IGRID . EQ . 3) THEN 
XD (1) = TH3 
XD (2) = CM3 
XD (3) = P3 
END IF 
END IF 
C 

DO 10 1=1, N 
C 

IF (XD (3) . NE .0.0. AND . XCI ( I ) ,LE.XD<3)) THEN 
C 

YC(I)=(XD(2) / (XD(3) *XD (3) ) ) * (2 . *XD (3) *XCI (I) -XCI (I) *XCI ( I) ) 
DYCDXCI(I) = (XD(2)/(XD(3)*XD(3)) ) * (2 . *XD (3) -2 . *XCI (I) ) 

C 

ELSEIF (XD (3) .NE. 0.0. AND. XCI (I) .GT .XD (3) ) THEN 
C 

YC (I) = (XD (2) / ( (1 . -XD (3) ) * (1 . -XD (3) ) ) ) * 

1 (1.-2. *XD (3) +2 . *XD (3) *XCI (I) -XCI (I) *XCI (I) ) 

DYCDXCI (I) = (XD (2) / ( (1 . -XD (3) ) * (1 . -XD (3))))*(2. *XD (3) -2 . *XCI (I) ) 

C 

ELSEIF (XD (3) .EQ. 0.0) THEN 


C 


YC (I) =0 . 0 


oo 



o o 


File: a Printed Tue Jan 01 00:00:00 198 Page: 7 


C 

10 


DYCDXCI (I) =0.0 

END IF 
CONTINUE 


49 


c 

c 


c 


c 

c 


c 

15 

C 

C 

C 

c 

c 


c 

25 

c 


DO 15 1=1, N 

IF(XCI(I) .LE. CHORD) THEN 
YT ( I) = (XD ( 1) /0 .2) * 

1 (A1*SQRT (XCI (I) ) +A2*XCI (I) +A3*XCI (I) *XCI (I) 

2 +A4*XCI (I) *XCI (I) *XCI (I) +A5*XCI (I) *XCI (I) *XCI (I) *XCI (I) ) 

DYTDXCI (I) = (XD (1) /0 .2) * (Al* (0 .5/SQRT (XCI (I) ) ) +A2+2 . *A3*XCI (I) 
1 +3 . *A4*XCI (I) *XCI (I) +4 . *A5*XCI (I) *XCI (I) *XCI (I) ) 

ELSE 


YT (I) =0 . 0 
DYTDXCI (I) =0.0 


END IF 
CONTINUE 


SURFACE COORDINATES (EQS . 5,6) 

DO 25 I = 1,N 

X(I) = XCI (I) 

IF ( II . EQ . 1) THEN 

Y ( I) = YC (I) + YT (I) 

DYDX(I) = DYCDXCI (I) + DYTDXCI (I) 
ANGLE ( I ) = ATAN (DYDX (I) ) 

ANGLE (1) = PI/2 

ELSE IF ( I I . EQ . 2 ) THEN 
Y (I) = YC (I) - YT (I) 

DYDX (I) = DYCDXCI (I) - DYTDXCI (I) 

ANGLE (I) = ATAN (DYDX (I) ) 

ANGLE (1) = PI/2 

ELSE 

WRITE (*, *) 9 Trouble in FOIL' 

STOP 
END IF 

CONTINUE 

RETURN 

END 



APPENDIX C 


FORTRAN LISTING OF NACA FOUR-DIGIT GRID 
SENSITIVITY : SENSIT 


50 
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SUBROUTINE DXYDXYD (XX, YY, NI , NJ, I CUR, FLAG, ISBNS, SRAN) 

PARAMETER ( ICHOICE=8 ) 

$ INCLUDE tbggl.inc 
$ INCLUDE tbggS.inc 

DIMENSION XX (NGDIM, NGDIM) , YY (NGDIM, NGDIM) 

COMMON/ XYDXYD/ DXDTP ' NGDIM, NGDIM) , DYDTH (NGDIM, NGDIM) , 

1 DXDOV JGDIM, NGDIM) , DYDCM (NGDIM, NGDIM) , 

2 DXDP (NGDIM, NGDIM) , DYDP (NGDIM, NGDIM) 

COMMON /DXYBDES/DXBD CM (NGDIM) , DYBDCM (NGDIM) , DXBDTH (NGDIM) , 

1 DYBDTH (NGDIM) , DXBDP (NGDIM) , DYBDP (NGDIM) , DXBEDTH (NGDIM) , 

2 DYBEDTH (NGDIM) , DXBEDCM (NGDIM) , DYBEDCM (NGDIM) , DXBEDP (NGDIM) , 

3 DYBEDP (NGDIM) 

COMMON/PQ/P (MGDIM) , Q (MGDIM) , T (MGDIM, MGDIM) 

COMMON/ JAYS/ J1 , J2 , J3 
LOGICAL FLAG (ICHOICE) 

C 

C 

C 

C 

C 

C GRID SENSITIVITY WITH RESPECT TO DESIGN VARIABLES 

C 

C WRITTEN BY : IDEEN SADREHAGHIGHI 

C 

Get Bottom Boundary Sensitivity 
C 

CALL DXYBDD(XX, YY,NI,NJ, ICUR, SRAN, CHEAT) 

C 

C Sensitivity Without Arc-length Distribution For Normal Direction 

C 

DO 10 J = 1 , NJ 
C 

DO 10 X « 1 , NX 
C 

CALL BLENDF (T (I, J) , ALFA1 , ALFA2 , ALFA3 , ALFA4 ) 

C 

C 

DXDTH (I, J) » ALFA1 * DXBDTH (I) + P(I) * ALFA3 * DXBEDTH (I) 

DYDTH (I, J) = ALFA1 * DYBDTH (I) + P(I) * ALFA3 * DYBEDTH (I) 

DXDCM(I,J) = ALFA1 * DXBDCM(I) + P(I) * ALFA3 * DXBEDCM ( I ) 

DYDCM (X, J) = ALFA1 * DYBDCM (I) + P(I) * ALFA3 * DYBEDCM ( I ) 

DXDP (I,J) = ALFA1 * DXBDP (I) + P(I) * ALFA3 * DXBEDP ( I) 

DYDP (I,J) = ALFA1 * DYBDP (I) + P(I) * ALFA3 * DYBEDP ( I) 

C 

10 CONTINUE 
C 

FLAG (ISBNS) = .TRUE. 

C 

RETURN 

END 

C 

c 

SUBROUTINE DXYBDD (XX, YY, NI, NJ, ICUR, SRAN, CHEAT) 

PARAMETER ( ICHOICE=8 ) 

$ INCLUDE tbggl.inc 

DIMENSION XX (NGDIM, NGDIM) , YY (NGDIM, NGDIM) , 

1 DXBUDTH (NGDIM) , DXBUDCM (NG^IM) , DXBUDP (NGDIM) , 

2 DYBUDTH (NGDIM) , DYBUDCM (NGDIM) , DYBUDP (NGDIM) , 

51 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


1 


c 


2 

c 

c 

c 

c 

c 

c 

c 


3 DXBLDTH (NGDXM) , DXBLDCM (NGDIM) , DXBLDP (NGDIM) , 

4 DYBLDTH (NGDIM) , DYBLDCM (NGDIM) , DYBLDP (NGDIM) , 

5 DXUEDTH (NGDIM) , DXUEDCM (NGDIM) , DXUEDP (NGDIM) , 

6 DYUEDTH (NGDIM) , DYUEDCM (NGDIM) , DYUEDP (NGDIM) , 

7 DXLEDTH (NGDIM) , D XX ED CM (NGDIM) , DXXEDP (NGDIM) , 

8 DYLEDTH (NGDIM) , DYLEDCM (NGDIM) , DYLEDP (NGDIM) , 

9 XU (NGDIM) , XL (NGDIM) , YU (NGDIM) , YL (NGDIM) 

COMMON /DXYBDES /DXBDCM (NGDIM) , DYBDCM (NGDIM) , DXBDTH (NGDIM) , 

1 DYBDTH (NGDIM) , DXBDP (NGDIM) , DYBDP (NGDIM) , DXBEDTH (NGDIM) , 

2 DYBEDTH (NGDIM) , DXBEDCM (NGDIM) , DYBEDCM (NGDIM) , DXBEDP (NGDIM) , 

3 DYBEDP (NGDIM) 

COMMON/ JAYS/ Jl, J2 , J3 

COMMON/TANGENT/THETAU (NGDIM) ; THETAL (NGDIM) , DYUDX (NGDIM) , DYLDX (NGDIM) , NFOIL 

COMMON / A/Al , A2 , A3 , A4 , A5 

COMMON /DESIGN/CM, P, TH, R, WLENGTH, CHORD 

COMMON /TEMP1/XNEW (NGDIM) , YNEW (NGDIM) , KNEW (NGDIM) 

LOGICAL CHEAT 


COMPUTES ANALITICALLY THE DERIVATIVE OF BOUNDARY COORDINATES 
WITH RESPECT TO DESIGN VARIABLES DXB/DXD 

«< DESIGN VARIBLES »> 


CM CAMBER 

P LOCATION OF CAMBER 

TH MAX. THICKNESS 


I UPPER = J2 

I LOWER = NI - J2 + 1 

IF ( I UPPER . NE . I LOWER) THEN 

WRITE (*,*) f ERROR FROM DXYBDD' 
STOP 
END IF 

DO 1 I « 1 , I UP PER 

XU(I) = XX ( J2-I+1, 1) - R 
YTJ(I) = YY ( J2-I+1, 1) 

CONTINUE 

DO 2 I = 1 , I LOWER 

XL (I) = XX ( I+J2-1 , 1) - R 
YL (I) = YY ( I+J2-l f 1) 

CONTINUE 


BOUNDARY DESIGN SENSITIVITY DERIVATIVES OF AIRFOIL 


52 
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C 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 


IDUMMY = J2 - J1 
DO 16 I = 1, IUPPER 

IF ( X . LE . IDUMMY) THEN 

CAMBERED AIRFOIL 

FORWARD OF MAX. CAMBER 

IF (P . NE .0.0 .AND .XU (I) . LE . P ) THEN 

YC = (CM/(P*P))*(2.*P*XU(I)-XU(I)*XU(I)) 

DYCDCM = (l./(P*P) )*(2.*P*XU(I)-XU(I)*XU(I)) 

DYCDP = (2 . *CM/ (P*P) ) * (-XU (I) +XU (I) *XU (I) /P) 

D2YCDCM= (1 . / (P*P) ) * (2 . *P-2 . *XU (I) ) 

D2YCDP = (2.*CM/(P*P))*(2.*XU(I)/P-1.) 

ELSE IF (P . NE .0.0. AND . XU (I) . GT . P ) THEN 

AFT OF MAX. CAMBER 

YC = (CM/ ( (1 . -P) * (1 . -P) ) ) * (1 . -2 . *P+2 . *P*XU (I) - 
XU(I) *XU(I) ) 

DYCDCM = (l./((l.-P)Ml.-P)))*(l-2.*P+2.*P*XU(I) 

-XU (I) *XU(I) ) 

FACT0R1 = CM/((1.-P)*(1.-P)) 

FACT0R2 = 2 . / (1 . -P) 

FACT0R3 = 1.-2. *P+2 . *P*XU (I) -XU (I) *XU (I) 

DYCDP = FACT0R1 * (FACT0R2 * FACT0R3 - 2. + 2.*XU(I)) 
D2YCDCM= (1 . / ( (1 . -P) * (1 . -P) ) ) * (2 . *P-2 . *XU(I) ) 

D2YCDP = (2 . *CM/ ((l.-P)*(l.-P)))*(2.*( (P-XU (I))/(l.-P))+l.) 

ELSEIF (CM . EQ .0.0. AND . P . EQ .0.0) THEN 

SYMMETRICAL AIRFOIL 

YC =0.0 
DYCDCM =0.0 
DYCDP =0.0 
D2YCDCM= 0.0 
D2YCDP =0.0 


END IF 


D2YTDTH = (1./0.2) *(0.5*A1/SQRT(XU(I))+A2 
+2 . *A3*XU(I) +3 . *A4*XU (I) *XU(I) +4 . *A5*XU(I) *XU(I) *XU (I) ) 

FACTOR = A1*SQRT (XU (I) ) +A2*XU (I) +A3*XU(I) *XU (I) 

+A4*XU (I) *XU(I) *XU ( I ) +A5 *XU ( I ) *XU(I) *XU(I) *XU(I) 

YT = (TH/0 .2) *FACT0R 

DYDTH = YT/TH 
DYBUDTH (I) = DYDTH 
DYBUDCM(I) = DYCDCM 
DYBUDP(I) = DYCDP 
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C 

C 


c 


c 

c 

c 


c 


c 

c 

c 


c 


c 

c 

c 


c 


TIMESU = COS (THETAU (I))/(l.+ (DYUDX (I) *DYUDX ( I) ) ) 

DXUEDTH ( I ) = + D2YTDTH * - TIMESU 


TIMESU = SIN (THETAU ( I) ) / ( 1 . + (DYUDX (I) *DYUDX ( I) ) ) 

DYUEDTH (I) = + D2YTDTH * - TIMESU 


EVALUATING DXETA/DCM , DYETA/DCM 


TIMESU = COS (THETAU (I) ) / (1 .+ (DYUDX (I) *DYUDX(I) ) ) 

DXUEDCM ( I ) = + D2YCDCM * - TIMESU 


TIMESU = S IN ( THETAU ( I ) ) / ( 1 . + (DYUDX ( I ) *DYUDX ( I ) ) ) 

DYUEDCM ( I ) = + D2YCDCM * - TIMESU 

. . EVALUATING DXETA/DP , DYETA/DP 


TIMESU = COS (THETAU (I) ) / (1 . + (DYUDX (I) *DYUDX ( I) ) ) 

DXUEDP (I) = + D2YCDP * - TIMESU 


TIMESU = SIN (THETAU (I) ) / (1 . + (DYUDX (I) *DYUDX (I) ) ) 

DYUEDP (I) = + D2YCDP * - TIMESU 


Singularity At Nose .... Slope dy/dx = Infinite 


IF ( I . EQ . 1 ) THEN 
DXUEDTH (I) =0.0 
DYUEDTH (I) =0.0 
DXUEDCM (I) = 0.0 
DYUEDCM (I) =0.0 
DXUEDP (I) =0.0 
DYUEDP (I) =0.0 
END IF 


ELSE 


C 

C SET BOUNDARY SENSITIVITY DERIVATIVES TO ZERO IN WAKE REGION 

C 

DXBUDCM ( I ) =0.0 
DYBUDCM (I) =0.0 
DXBUDP (I) =0.0 
DYBUDP (I) =0.0 
DXBUDTH (I) =0.0 
DYBUDTH ( I ) =0.0 
DXUEDTH (I) =0.0 
DYUEDTH (I) =0.0 
DXUEDCM ( I ) = 0.0 
DYUEDCM (I) =0.0 
DXUEDP (I) =0.0 
DYUEDP (I) = 0.0 
C 


C 

16 

C 

C 


END IF 
CONTINUE 

IL = J2 - J1 + 1 


C 

C 
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DXB/DXD = DXB/DR * DR/DXD 


C 


C 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 


CALL DXBDXD ( 1 , IL, XU, YU, DXBUDCM, DXBUDTH, DXBUDP, CHEAT) 

IDUMMY = J3 - J2 
DO 17 I = 1, I LOWER 

IF ( I . LE . IDUMMY) THEN 

CAMBERED AIRFOIL 

FORWARD OF MAX. CAMBER 

IF (P .NE .0.0. AND . XL (I) . LE . P) THEN 

YC = (CM/ (P*P) ) * (2 . *P*XL (I) -XL (I) *XL (I) ) 

DYCDCM = (l./(P*P))*(2. *P‘XL (I) -XL (I) *XL (I) ) 

DYCDP ■ (2.*CM/(P*P))*(-XL(I)+XL(I)*XL(I)/P) 

D2YCDCM= (l./(P*P))*(2. *P-2 . *XL (I) ) 

D2YCDP = (2 . *CM/ (P*P) ) * (2 . *XL (I) /P-1 . ) 

ELSE IF (P . NE .0.0. AND . XL (I) . GT . P ) THEN 

AFT OF MAX. CAMBER 

YC = (CM/ ( (1 . -P) * (1 . -P) ) ) * (1 . -2 . *P+2 . *P*XL(I) - 
XL (I) ‘XL (I) ) 

DYCDCM = (1 . / ( (1 . -P) * (1 . -P) ) ) * (1 . -2. ‘P+2 . *P*XL(I) 

-XL (I) *XL (I) ) 

FACTOR1 = CM/ ( (1. -P) * (1. -P) ) 

FACTOR2 = 2 . / (1 . -P) 

FACTOR3 = 1 . -2 . *P+2 . *P‘XL (I) -XL (I) ‘XL (I) 

DYCDP * FACTOR1 * (FACTOR2 * FACTOR3 - 2. + 2. ‘XL (I)) 
D2YCDCM= (1 . / ( (1 . -P) * (1 . -P) ) ) * (2 . ‘P-2 . ‘XL (I) ) 

D2YCDP = (2 . *CM/ ( (1 . -P) * (1 . -P) ) ) * (2 . * ( (P-XL(I) ) / (1 . -P) ) +1 . ) 

ELSEIF (CM . EQ .0.0. AND . P . EQ .0.0) THEN 

SYMMETRICAL AIRFOIL 

YC =0.0 

DYCDCM =0.0 
DYCDP =0.0 
D2YCDCM= 0.0 
D2YCDP =0.0 


C 

c 


c 


c 

c 


END IF 


D2YTDTH = (1 . /0 . 2) * (0 . 5*A1/SQRT (XL (I) ) +A2 

+2 . *A3*XL (I) +3 . *A4*XL (I) *XL (I) +4 . *A5*XL (I) *XL (I) *XL(I) ) 

FACTOR = A1*SQRT (XL (I) ) +A2*XL (I) +A3*XL (I) *XL ( I) 

+A4*XL(I) *XL (I) *XL (I) +A5*XL (I) *XL(I) *XL(I) *XL(I) 

YT = (TH/0 .2) ‘FACTOR 

DYDTH = YT/TH 
DYBLDTH (I) = - DYDTH 
DYBLDCM ( I ) = DYCDCM 
DYBLDP(I) = DYCDP 
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C 

c 

c 

c 


c 


c 

c 

c 


c 


c 

c 

c 


c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

17 

c 

c 

c 


EVALUATING DXETA/DTH , DYETA/DTH 


TIMESL = COS ( THETAL ( I ) ) / ( 1 . + (DYLDX ( I ) *DYLDX ( I ) ) ) 

DXLEDTH (I) = - D2YTDTH * + TIMESL 


TIMESL = S IN ( THETAL ( I ) ) / ( 1 . + (DYLDX ( I ) *DYLDX ( I ) ) ) 

DYLEDTH (I) = - D2YTDTH * + TIMESL 


EVALUATING DXETA/DCM , DYETA/DCM 


TIMESL = COS ( THETAL ( I ) ) / ( 1 . + (DYLDX ( I ) *DYLDX ( I ) ) ) 

DXLEDCM ( I ) = 4- D2YCDCM * + TIMESL 


TIMESL = SIN (THETAL (I))/{1.+ (DYLDX (I) *DYLDX (I) ) ) 

DYLEDCM (I) = + D2YCDCM * + TIMESL 


EVALUATING DXETA/DP , DYETA/DP 


TIMESL = COS { THETAL ( I ) ) / ( 1 . + (DYLDX ( I ) *DYLDX ( I ) ) ) 

DXLEDP (I) = + D2YCDP * + TIMESL 


TIMESL = SIN (THETAL ( I) ) / (1 . + (DYLDX (I) *DYLDX ( I) ) ) 

DYLEDP(I) = + D2YCDP * + TIMESL 


Singularity At Nose .... Slope dy/dx ~ Infinite 


IF ( I . EQ . 1) THEN 
DXLEDTH (I) =0.0 
DYLEDTH (I) =0.0 
DXLEDCM (I) =0.0 
DYLEDCM (I) =0.0 
DXLEDP (I) =0.0 
DYLEDP (I) = 0.0 
END IF 


ELSE 

SET SENSITIVITY DERIVATIVES TO ZERO IN WAKE REGION 


DXBLDCM ( I ) 

= 

0.0 

DYBLDCM ( I ) 

= 

0.0 

DXBLDP (I) 

= 

0.0 

DYBLDP ( I ) 

= 

0.0 

DXBLDTH (I) 

= 

0.0 

DYBLDTH (I) 

= 

0.0 

DXLEDTH (I) 

= 

0.0 

DYLEDTH (I) 

= 

0.0 

DXLEDCM (I) 

= 

0.0 

DYLEDCM (I) 

=5 

0.0 

DXLEDP (I) 

= 

0.0 

DYLEDP (I) 

= 

0.0 

END IF 



CONTINUE 
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C 

C 

c 


c 

50 

C 

C 


c 

60 

C 


C 


CALL DXBDXD (2, IL, XL, YL, DXBLDCM, DXBLDTH, DXBLDP, CHEAT) 


DO 50 I = 1 , IUPPER 

DXBDTH (IUPPER - I + 1) 
DYBDTH( IUPPER - I + 1) 
DXBDCH ( IUPPER - I + 1) 
DYED CM ( IUPPER - I + 1) 
DXBDP (IUPPER - I + 1) 
DYBDP (IUPPER - I + 1) 
DXBEDTH (IUPPER- I + 1) 
DYBEDTH( IUPPER- I + 1) 
DXBEDCM( IUPPER- I + 1) 
DYBEDCM( IUPPER- I + 1) 
DXBEDP (IUPPER - I + 1) 
DYBEDP ( IUPPER - I + 1) 

CONTINUE 

DO 60 I = 1, I LOWER 

DXBDTH (IUPPER - 1 + I) 

DYBDTH( IUPPER - 1 + I) 

DXBDCM ( IUPPER - 1 + I) 

DYBDCM( IUPPER - 1 + I) 

DXBDP (IUPPER - 1 + I) 

DYBDP (IUPPER - 1 + I) 

DXBEDTH (IUPPER- 1 + I) 
DYBEDTH (IUPPER- 1 + I) 
D XB ED CM ( IUPPER- 1 + I) 
DYBEDCM( IUPPER- 1 + I) 
DXBEDP (IUPPER - 1 + I) 

DYBEDP (IUPPER - 1 + I) 


DXBUDTH (I) 
DYBUDTH (I) 
DXBUDCM ( I ) 
DYBUDCM (I) 
DXBUDP (I) 
DYBUDP (I) 
DXUEDTH (I) 
DYUEDTH (I) 
DXUEDCM (I) 
DYUEDCM(I) 
DXUEDP (I) 
DYUEDP (I) 


DXBLDTH (I) 
DYBLDTH (I) 
DXBLDCM (I) 
DYBLDCM (I) 
DXBLDP (I) 
DYBLDP (I) 
DXLEDTH (I) 
DYLEDTH (I) 
DXLEDCM (I) 
DYLEDCM(I) 
DXLEDP (I) 
DYLEDP (I) 


CONTINUE 


RETURN 

END 


C** ********************************************************* *********** 


C 

SUBROUTINE DXBDXD(II, I L, XNEW, YNEW, DXBDCM, DXBDTH, DXBDP, CHEAT) 

$ INCLUDE tbggl.inc 

DIMENSION XNEH(NGDIM) , YNEW (NGDIM) , DXBDCM (NGDIM) , DXBDTH (NGDIM) 

1 , DXBDP (NGDIM) , ROLD (NGDIM) , RNEW (NGDIM) , XOLD (NGDIM) 

2 , YOLD (NGDIM) , DYDTH (NGDIM) , DYDCM (NGDIM) , DYDP (NGDIM) 

3 , DRDCM (NGDIM) , DRDTH (NGDIM) , DRDP (NGDIM) , DRDCMN (NGDIM) 

4 , DRDTHN (NGDIM) , DRDPN (NGDIM) 

COMMON/OLD/XUS (NGDIM) , YUS (NGDIM) , XLS (NGDIM) , YLS (NGDIM) 

COMMON/TANGENT/THETAU (NGDIM) , THETAL (NGDIM) , DYUDX (NGDIM) , DYLDX (NGDIM) ,NFOIL 

COMMON/ A/ A1 , A2 , A3 , A4 , A5 

COMMON /DESIGN/CM, P, TH, R, WLENGTH, CHORD 

COMMON /DELTAXD/TH1, CM1, PI, TH3, CM3, P3 

COMMON/DUMMY/X (NGDIM) , Y (NGDIM) , XI (NGDIM) , X3 (NGDIM) , RO (NGDIM) , RI (NGDIM) , 

1 XDIST (NGDIM) , DUM1 (NGDIM) , DUM2 (NGDIM) , XCI (NGDIM) , ROBAR (NGDIM) 

LOGICAL CHEAT 
DATA CK, DELTA/4. 0,0. 0001/ 


C 

C 
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C 


C 

c 


c 

c 



c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


goto 999 

DO 10 I = 1 , NGDIM 

IF (II. EQ . 1 ) THEN 

XOLD (I) =XUS (I) 
YOLD (I) =YUS (I) 

ELSE 

XOLD ( I ) =XLS ( I ) 
YOLD ( I ) =YLS ( I ) 

END IF 

CONTINUE 


DO 15 I = 1, NGDIM 
IF (XOLD (I) .LE. CHORD) THEN 

CAMBERED AIRFOIL 

FORWARD OF MAX. CAMBER 

IF (P.NE. 0.0. AND. XOLD (I) . LE . P ) THEN 

YC ■ (CM/(P*P) )*(2.*P*XOLD(I)-XOLD(I)*XOLD(I)) 

DYCDCM = (1 . / (P*P) ) * (2 . *P*XOLD(I) -XOLD (I) *XOLD(I) ) 
DYCDP a (2 . *CM/ (P*P) ) * (-XOLD (I) +XOLD(I) *XOLD(I) /P) 

ELSEIF (P . NE . 0 . 0 . AND . XOLD (I) . GT . P ) THEN 

AFT OF MAX. CAMBER 

YC = (CM/ ( (l.-P) * (l.-P) ) ) *(1.-2 . *P+2.*P*XOLD(I) - 
XOLD (I) *XOLD (I) ) 

DYCDCM = (l./((l.-P)*(l.-P) ) ) * (1 . -2 . *P+2 . *P*XOLD (I) 

-XOLD (I) *XOLD (I) ) 

FACTOR1 = CM/((1.-P)*(1.-P)) 

FACTOR2 » 2 . / (1 . -P) 

FACTOR3 a 1 . -2 . *P+2 . *P*XOLD (I) -XOLD (I) *XOLD (I) 

DYCDP = FACTOR1 * (FACTOR2 * FACTOR3 - 2. + 2.*XOLD(I)) 

ELSEIF (CM . EQ . 0 . 0 . AND . P . EQ . 0 . 0 ) THEN 

SYMMETRICAL AIRFOIL 

YC =0.0 
DYCDCM =0.0 
DYCDP =0.0 


END IF 
ELSE 

WRITE ( * , * ) ' ERROR FROM DXBDXIgg XOLD IS OUT OF DOMAIN' 

STOP 


o o ooo non on ooo ooo no 
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END IF 


FACTOR = Al*SQRT(XOLD(I) ) +A2*XOLD (I) +A3*XOLD (I) *XOLD (I) 
+A4*XOLD (I) *XOLD (I) *XOLD (I) +A5*XOLD(I) *XOLD (I) *XOLD (I) *XOLD (I) 

YT = (TH/0 .2) *FACTOR 

IF(II.EQ.2)YT=-YT 

DYDTH (I) = YT/TH 
DYDCM(I) = DYCDCM 
DYDP(I) = DYCDP 

15 CONTINUE 


CALL ARC (XOLD, YOLD, ROLD, NGDIM, RMAX) 
CALL ARC (XNEW, YNEW, RNEW, IL , RMAX) 


DRDCM(l) =0 . 0 
DRDTH (1) =0 . 0 
DRDP (1) =0 . 0 

DO 20 I = 2 , NGDIM 

FACTOR = (YOLD (I) - YOLD (1-1) )/ (ROLD (I) -ROLD (1-1) ) 

DRDCM(I) = FACTOR * (DYDCM(I) -DYDCM(I-l) ) + DRDCM(I-l) 

DRDTH (I) = FACTOR * (DYDTH (I) -DYDTH (1-1) ) + DRDTH (1-1) 

DRDP (I) = FACTOR * (DYDP (I)-DYDP (1-1) ) + DRDP (1-1) 

20 CONTINUE 


CALL INTERPO (ROLD, DRDCM, NGDIM, RNEW, DRDCMN, IL) 
CALL INTERPO (ROLD, DRDTH, NGDIM, RNEW, DRDTHN, IL) 
CALL INTERPO (ROLD, DRDP , NGDIM, RNEW, DRDPN , IL) 


C 

C 

C 

C 


DO 25 I = 1 , IL 
IF(I.EQ.l) THEN 

DXBDR a (XNEW (1+1) -XNEW (I) ) / (RNEW (1+1) -RNEW (I) ) 
ELSEIF ( I . EQ . IL) THEN 

DXBDR = (XNEW(I) -XNEW(I-l) ) / (RNEW(I) -RNEW(I-l) ) 

ELSE 

DXBDR = (XNEW ( 1+1) -XNEW ( I — 1 ) ) / (RNEW (1+1) -RNEW (I- 1) ) 
END IF 

IF ( I I . EQ . 1 ) DYDX=DYUDX ( I ) 

IF ( II . EQ . 2 ) DYDX=DYLDX ( I ) 

FACTOR=SQRT(l. + DYDX*DYDX) 

DXBDR=1 . /FACTOR 

59 

DXBDCM(I) = DXBDR * DRDCMN (I) 


o o 
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DXBDTH(I) = DXBDR * DRDTHN(I) 

DXBDP (I) = DXBDR * DRDPN (I) 

C 

25 CONTINUE 
C 

999 CONTINUE 
goto 1099 
IF (CHEAT) THEN 
C 

TH1 = TH + DELTA 
TH3 = TH - DELTA 
IF (CM. NE. 0.0) THEN 
CM1 = CM + DELTA 
CM3 = CM - DELTA 
ELSE 
CM1=0 . 0 
CM3=0 .0 
END IF 

IF (P .NE .0.0) THEN 
PI = P + DELTA 
P3 = P - DELTA 
ELSE 
P1=0 .0 
P3=0 . 0 
END IF 
C 

IGRID = 1 
1000 CONTINUE 
C 

DO 35 I— 1, NGDIM 
C 

C = FLOAT (1-1) /FLOAT (NGDIM-1) 

XCI(I) = ( (EXP (CK*C) -1 . ) / (EXP (CK) -1 . ) ) *CHORD 
C 

35 CONTINUE 

C 

CALL FOIL (XCI,X,Y,DUM1,DUM2, NGDIM, II, CHEAT, IGRID) 
CALL ARC (X,Y,RI, NGDIM, RMAX) 

CALL BATCH (IL, ROBAR, 1) 

C 

DO 45 1=1, IL 

RO ( I ) = ROBAR ( I ) *RMAX 
45 CONTINUE 

C 

CALL INTERPO (RI,X, NGDIM, RO,XDIST, IL) 

CALL FOIL (XDIST, X, Y, DUM1, DUM2 , IL, II, CHEAT, IGRID) 

C 

IF ( IGRID . EQ . 1 ) THEN 
DO 65 1=1, IL 
X1(I) = X (I) 

65 CONTINUE 

IGRID = 3 
GOTO 1000 

ELSEIF ( IGRID . EQ . 3 ) THEN 
DO 85 1=1, IL 
X3(I) = X (I) 

85 CONTINUE 

END IF 

60 
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Finite Differencing 


DELTATH 

DELTACM 

DELTAP 


TH1 

CM1 

PI 


- TH3 

- CM3 

- P3 


DO 75 1=1, IL 

DXBDTH(I) = (XI (I) - X3 ( I) ) / (TH1 - TH3) 
IF (DELTACM. NE. 0.0) THEN 

DXBDCM(I) = (XI (I) -X3 (I) ) / (CM1 - CM3) 

ELSE 

DXBDCM(I) = 0.0 
END IF 

IF (DELTAP . NE . 0 . 0 ) THEN 


75 


DXBDP(I) 
ELSE 
DXBDP(I) 
END IF 
CONTINUE 


(XI (I) - X3 (I) ) / (PI - P3) 


= 0.0 


END IF 


1099 continue 

do 2000 1=1, IL 
DXBDTH(I) = 
DXBDCM (I) = 
DXBDP (I) = 
2000 CONTINUE 
RETURN 
END 


0.0 

0.0 

0.0 
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